Manejando datos espaciales con R - Parte I

[🔴] 1. Instalación de paquetes necesarios

pkg <- c(
  "sf",        # Paquete para manejar datos espaciales
  "tidyverse", # Paquete para ciencia de datos
  "mapview",   # Paquete para visualizar de forma interactiva datos espaciales
  "raster",    # Paquete para manejar datos raster
  "lattice",   # Paquete para la visualización de datos 
  "cptcity",   # Paquete para manejar paleta de colores
  "tmap"   ,   # Paquete para la elaboración de mapas temáticos
  "pacman" ,   # Para activar varios paquetes a la misma vez
  "ggspatial", # Complementos adicionales para la elaboración de mapas
  "cowplot")   # Paquete para mesclar plots

install.packages(
 pkgs = pkg,
 dependencies = TRUE   
)

[🔴] 2. Activación de los paquetes necesarios

library(pacman)
p_load(
  sf, tidyverse, mapview, raster,
  lattice, cptcity, tmap, ggspatial,
  cowplot
  )

[🔴] 3. Lectura datos vectoriales

loreto <- st_read('datos/vector/distritos_loreto.gpkg') #read_sf()

Reading layer `distritos_loreto' from data source 
  `/home/barja/Escritorio/materiales/markdown/datos/vector/distritos_loreto.gpkg' 
  using driver `GPKG'
Simple feature collection with 53 features and 8 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -77.82596 ymin: -8.716139 xmax: -69.94904 ymax: -0.03860597
Geodetic CRS:  WGS 84

class(loreto) # Indetificar el tipo de objeto  

[1] "sf"         "data.frame"

head(loreto)  # Visualizar las 6 primeras filas 

Simple feature collection with 6 features and 8 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -76.63459 ymin: -8.716139 xmax: -72.85497 ymax: -4.970566
Geodetic CRS:  WGS 84
  IDDIST      NOMBDIST      NOMBPROV NOMBDEP      NOM_CAP   AREA_KM. Pob2017
1 160601     CONTAMANA       UCAYALI  LORETO    CONTAMANA 11029.0188   23883
2 160505        MAQUIA       REQUENA  LORETO SANTA ISABEL  4891.5766    7304
3 160603 PADRE MARQUEZ       UCAYALI  LORETO     TIRUNTAN  2434.9196    3697
4 160602      INAHUAYA       UCAYALI  LORETO     INAHUAYA   629.6961    1738
5 160205       JEBEROS ALTO AMAZONAS  LORETO      JEBEROS  4665.0924    3900
6 160511     YAQUERANA       REQUENA  LORETO      ANGAMOS 11049.2448    1929
  FUENTE                           geom
1   INEI MULTIPOLYGON (((-74.88129 -...
2   INEI MULTIPOLYGON (((-74.49483 -...
3   INEI MULTIPOLYGON (((-74.64314 -...
4   INEI MULTIPOLYGON (((-75.10987 -...
5   INEI MULTIPOLYGON (((-75.91967 -...
6   INEI MULTIPOLYGON (((-72.85497 -...

tail(loreto)  # Visualizar las últimas 6 filas

Simple feature collection with 6 features and 8 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -77.73971 ymin: -5.746232 xmax: -70.06006 ymax: -0.03860597
Geodetic CRS:  WGS 84
   IDDIST                NOMBDIST          NOMBPROV NOMBDEP
48 160110          TORRES CAUSANA            MAYNAS  LORETO
49 160803 TENIENTE MANUEL CLAVERO          PUTUMAYO  LORETO
50 160702              CAHUAPANAS DATEM DEL MARAÑON  LORETO
51 160701                BARRANCA DATEM DEL MARAÑON  LORETO
52 160802            ROSA PANDURO          PUTUMAYO  LORETO
53 160804                  YAGUAS          PUTUMAYO  LORETO
                     NOM_CAP  AREA_KM. Pob2017 FUENTE
48                   PANTOJA  6926.423    4230   INEI
49             S0PLIN VARGAS  9685.187    2317   INEI
50 SANTA MARIA DE CAHUAPANAS  4781.949    6336   INEI
51               SAN LORENZO  7372.721   12742   INEI
52            SANTA MERCEDES  7064.397     520   INEI
53                   REMANSO 18073.543    1277   INEI
                             geom
48 MULTIPOLYGON (((-74.8737 -0...
49 MULTIPOLYGON (((-75.10004 -...
50 MULTIPOLYGON (((-76.45404 -...
51 MULTIPOLYGON (((-75.91967 -...
52 MULTIPOLYGON (((-74.80014 -...
53 MULTIPOLYGON (((-70.93735 -...

names(loreto) # Identificar los nombres de las columnas

[1] "IDDIST"   "NOMBDIST" "NOMBPROV" "NOMBDEP"  "NOM_CAP"  "AREA_KM." "Pob2017" 
[8] "FUENTE"   "geom"    

dim(loreto)   # Total de filas x columnas

[1] 53  9

st_crs(loreto)[[1]] # SRC

[1] "WGS 84"

[🔴] 4. Visualización de vectoriales

# Usando las funciones nativas de R
plot(
  loreto["Pob2017"],
  axes = TRUE,
  border = 'grey'
  ) 

# Usando la sintaxis de tmap
loreto %>% 
  tm_shape() +
  tm_polygons(
    col = "Pob2017",
    palette = cpt(pal = "mpl_viridis"),
    border.col = "white",
    lwd = 0.5
    )

# Usando la sintaxis de ggplot
loreto %>% 
  ggplot() + 
  geom_sf(aes(fill = Pob2017)) + 
  scale_fill_gradientn(
    colours = cpt(pal = "mpl_viridis")
    ) + 
  theme_minimal()

# Visualización interactiva
loreto %>% 
  mapview(z = "AREA_KM.")

[🔴] 5. Geoprocesamiento con datos vectoriales

# Lectura de datos csv 
covid_data <- read_csv('datos/csv/data_positivos_covid.csv')
covid_positivo <- covid_data %>% 
  group_by(UBIGEO) %>% 
  summarise(casos_positos = n())

# Verificar que campos no se unen 
anti_join(
  loreto,
  covid_positivo,
  by = c("IDDIST"="UBIGEO")
)

Simple feature collection with 1 feature and 8 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -74.19259 ymin: -6.484717 xmax: -73.83009 ymax: -5.574155
Geodetic CRS:  WGS 84
  IDDIST NOMBDIST NOMBPROV NOMBDEP NOM_CAP AREA_KM. Pob2017 FUENTE
1 160509  TAPICHE  REQUENA  LORETO  IBERIA 2090.124     881   INEI
                            geom
1 MULTIPOLYGON (((-73.95703 -...

# Join por campo en común 
loreto_covid <- left_join(
  loreto,
  covid_positivo,
  by = c("IDDIST"="UBIGEO")
)

names(loreto_covid)

 [1] "IDDIST"        "NOMBDIST"      "NOMBPROV"      "NOMBDEP"      
 [5] "NOM_CAP"       "AREA_KM."      "Pob2017"       "FUENTE"       
 [9] "casos_positos" "geom"         

class(loreto_covid)

[1] "sf"         "data.frame"

# Generar areas de influencia 
# Reproyectar -> Centroides -> Buffer -> Divolse -> Visualizar resultado

loreto_centroid <- loreto_covid %>% 
  filter(casos_positos >= 100) %>% 
  st_transform(crs = 32718) %>% 
  st_centroid() %>% 
  st_buffer(dist = 50000) %>% 
  group_by() %>% 
  summarise() %>% 
  mutate(zona_influencia = unlist(st_area(geom))/1000000)

loreto_centroid %>% 
  mapview()

# Intersectar zonas de influencias con los distritos 

zona_riesgo <- st_intersection(
  x = loreto_centroid,
  y = loreto %>% st_transform(32718)
  )
zona_riesgo %>% 
  tm_shape() +
  tm_polygons("Pob2017")

# Disolver capas distritos -> provincias | distritos -> departameto
m1 <- loreto %>% mapview(z="Pob2017",layers.control.pos = "bottom")
provincias_loreto <- loreto %>% 
  group_by(NOMBPROV) %>% 
  summarise(Pob2017 = sum(Pob2017))

m2 <- provincias_loreto %>%
  mapview(z="Pob2017")

m1 | m2

[🔴] 6. Escritura de datos vectoriales

write_sf(
  zona_riesgo,
  "datos/output/zona_riesgo.gpkg",
  delete_dsn = TRUE
  )

[🔴] 7. Lectura con datos raster

# En los datos en raster existe tres tipos de datos: 
# 1. Layer: Raster con una sola bandas
# 2. Stack: Raster con bandas separadas
# 3. Brick: Raster con multiple bandas

# Lectura de un modelo de elevación 
dem <- raster('datos/raster/dem/elevacion.tif')
class(dem)

[1] "RasterLayer"
attr(,"package")
[1] "raster"

crs(dem)

CRS arguments: +proj=longlat +datum=WGS84 +no_defs 

res(dem)

[1] 0.0002777778 0.0002777778

extent(dem)

class      : Extent 
xmin       : -75.96014 
xmax       : -75.78153 
ymin       : -11.76347 
ymax       : -11.60264 

maxValue(dem)

[1] 32767

minValue(dem>=0)

[1] 1

summary(dem, maxsamp = ncell(dem))

        elevacion
Min.         3569
1st Qu.      4066
Median       4304
3rd Qu.      4491
Max.         4753
NA's            0

# Lectura de imagenes satelitales por multiples bandas
escenas <- list.files(
  path = 'datos/raster/bandas/',
  pattern = "*.tif"
  )
iquitos_stack <- stack(
  sprintf(
    "datos/raster/bandas/%s",
    escenas
    )
  )
class(iquitos_stack)

[1] "RasterStack"
attr(,"package")
[1] "raster"

iquitos_stack

class      : RasterStack 
dimensions : 1503, 1385, 2081655, 4  (nrow, ncol, ncell, nlayers)
resolution : 8.983394e-05, 8.983367e-05  (x, y)
extent     : -73.31241, -73.18799, -3.79745, -3.66243  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs 
names      :       b1,       b2,       b3,       b4 
min values : 186.5000, 315.0000, 233.1667, 284.0000 
max values :     6040,     5776,     5428,     5843 

# Lectura de satelitales como uno solo

iquitos_brick <- brick("datos/raster/sentinel2/Iquitos.tif")
class(iquitos_brick)

[1] "RasterBrick"
attr(,"package")
[1] "raster"

iquitos_brick

class      : RasterBrick 
dimensions : 1503, 1385, 2081655, 4  (nrow, ncol, ncell, nlayers)
resolution : 8.983153e-05, 8.983153e-05  (x, y)
extent     : -73.31241, -73.18799, -3.797448, -3.662431  (xmin, xmax, ymin, ymax)
crs        : +proj=longlat +datum=WGS84 +no_defs 
source     : Iquitos.tif 
names      : B4, B3, B2, B8 

[🔴] 8. Visualización de datos raster

# Funciones nativas de R
dem %>% plot()

# Funciones con Tmap
p2 <- dem %>% tm_shape() +
  tm_layout(legend.outside = T) + 
  tm_raster(palette = cpt(pal = "grass_elevation"))
p2

# Funciones con ggplot
p3 <- dem %>% as.data.frame(xy=TRUE) %>% 
  ggplot(aes(x = x,y = y ,fill = elevacion)) + 
  geom_tile() + 
  scale_fill_gradientn(colours = cpt(pal = "grass_elevation")) + 
  theme_minimal()
p3

# Visualización en RGB de Iquitos
iquitos_brick %>% 
  plotRGB(4,3,2,stretch = "lin")

# Visualización interactiva

iquitos_brick %>% 
  mapview::viewRGB(r = 4,g =3,b = 2 )

[🔴] 9. Geoprocesamientos con datos raster

# Corte por mascara 

iquitos_shp <- loreto %>% 
  filter(NOMBDIST == "IQUITOS")
corte <- iquitos_brick %>% 
  crop(iquitos_shp) %>% 
  mask(iquitos_shp)
plot(corte)

plotRGB(corte,4,3,2,stretch = "lin")

Operaciones algebráicas con bandas espectrales

• Objetivo: Analizar las respuestas espectrales de las diferentes

tipos de cobertura mediante las imágenes satelitales.

• NDVI : Normalized Difference Vegetation Index, índice que se utiliza para estimar la cantidad,calidad y desarrollode la vegetación.

• SAVI : Soil Adjusted Vegetation Index, una mejora del ndvi para evitar las distorciones en cuanto al análisis de vegetación que se encuentra sobre suelos expuestos.

Mas información click aquí

Entonces para calular en NDVI necesitamos la banda roja y la banda infraroja, en las imágenes sentinel2 estas son: * Banda roja : Banda 4 * Banda infraroja: Banda 8

ndvi <- function(nir,red){
  ndvi <- (nir - red)/(nir + red)
  return(ndvi)
}

iquitos_ndvi <- ndvi(
  nir = iquitos_brick$B8,
  red = iquitos_brick$B4
  )
pal <- cpt(pal = "grass_ndvi",rev = FALSE)
mapview(
  iquitos_ndvi,
  col.regions = pal,
  legend = TRUE
  )

Finalmente, para calular en SAVI necesitamos la banda roja y la banda infraroja, en las imágenes sentinel2 estas son: * Banda roja : Banda 4 * Banda infraroja: Banda 8 * factor L: 0.5

savi <- function(nir,red){
  savi <- ((nir - red)/(nir + red + L))*(1+L)
  return(savi)
}

iquitos_savi <- ndvi(
  nir = iquitos_brick$B8,
  red = iquitos_brick$B4
  )
mapview(
  iquitos_savi,
  legend = TRUE
  )

[🔴] 10. Escritura de datos raster

writeRaster(
  corte,
  "datos/output/iquitos_corte.tif",
  overwrite=TRUE
  )
writeRaster(
  iquitos_ndvi,
  "datos/output/iquitos_ndvi.tif",
  overwrite=TRUE
  )